R and geodata, part 2

Andreja Radović, PhD

25th of January 2018

Gridded data in ‘sp’

Again, we can create objects,here start from grid topology definition

[1] “SpatialGrid” attr(,“package”) [1] “sp”

Formal class ‘SpatialGrid’ [package “sp”] with 3 slots ..@ grid :Formal class ‘GridTopology’ [package “sp”] with 3 slots .. .. ..@ cellcentre.offset: num [1:3] 1 1 2 .. .. ..@ cellsize : num [1:3] 1 1 1 .. .. ..@ cells.dim : int [1:3] 3 4 6 ..@ bbox : num [1:3, 1:2] 0,5 0,5 1,5 3,5 4,5 7,5 .. ..- attr(*, “dimnames”)=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] “min” “max” ..@ proj4string:Formal class ‘CRS’ [package “sp”] with 1 slot .. .. ..@ projargs: chr NA

What we have made so far

Create Grid object from points

Start to build objects from grid topology, coordinates or read gridded data

   x y
1  1 1
2  2 1
3  3 1
4  1 2
5  2 2
6  3 2
7  1 3
8  2 3
9  3 3
10 1 4
11 2 4
12 3 4
[1] "SpatialPoints"
attr(,"package")
[1] "sp"

Visualise prepared grid

Continuation to full grid

[1] "SpatialPoints"
attr(,"package")
[1] "sp"
[1] FALSE

Continuation to full grid

[1] 1 2 3
[1] "data.frame"
[1] "SpatialGridDataFrame"
attr(,"package")
[1] "sp"

visualise prepared grid

Visaualisation with lattice graphics

Import gridded data with sp package

Common way of importing gridded data; calculate summary statistics…

dem_1k.asc has GDAL driver FITS and has 21 rows and 21 columns

Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 

0,000 6,419 157,996 301,073 511,735 1255,097

With all types of spatial data, statistical procedures from R

As example plot histogram of elevation values

Visualise gridded data with attributes

Export as …

ESRI asc, SAGA sgrd…

Gridded data with raster package

  • Extremely popular for analysis of big rasters

  • Data does not need to be in memory

  • Easy shift to ‘sp’ classes

  • Extract wanted spatial extent or aggregate by polygon, select on points…

Get data from Worldclim database - raster object

Plot raster object

Extract values for polygons

Values of Digital elevation model for each county

dem_1k.asc has GDAL driver FITS and has 21 rows and 21 columns

[[1]] [1] 722,6663 261,4872 370,5971 773,7128 121,7043 307,7529 601,9179 259,0665

[[2]] [1] 157,0298 122,3399 136,0429 153,9509 157,9960 149,9645 356,1018

[[3]] [1] 125,55239 143,63921 144,06564 97,78848

[[4]] [1] 98,8438 387,1294

[[5]] [1] 196,9077 168,1451 122,0084

[[6]] [1] 185,08521 377,42798 163,93018 281,47314 64,61814 173,73517 24,62883

[[7]] [1] 145,2687 143,5296 336,7780 201,0698 206,9406 757,8558 546,8431 [8] 352,1976 266,6391 691,3894

[[8]] [1] 198,0352 150,9576 123,2260 138,8395 153,0450

[[9]] [1] 285,3583 284,5465 202,5881 237,4919

[[10]] [1] 524,9102 652,7065 429,8729 710,6260 729,8898 779,5086 760,5005 [8] 746,7065 675,1872 809,3881 826,6930 702,2512 725,9460

[[11]] NULL

[[12]] [1] 92,74722 90,60442 84,34888 80,32372 142,12306 95,99942 78,94941 [8] 115,79662

[[13]] [1] 176,4084 382,2073 293,2750 217,3138 196,9911

[[14]] [1] 964,98822 497,56903 457,62909 807,84265 762,32831 171,38272 386,63550 [8] 97,36595 48,06548 26,26266

[[15]] [1] 95,6325 126,6646 141,5341 131,3787 105,1319 111,2801 188,6825 [8] 269,6385 188,2933 125,9783 133,8734 272,3000

[[16]] [1] 714,4169 371,5448 444,0142 626,7087 280,6295 547,9183 355,3295

[[17]]

228,8269

[[18]] [1] 114,6883 191,4677 153,4848 103,0313 412,7148 267,1815

[[19]] [1] 81,57622 75,48589

[[20]] [1] 831,36621 800,62671 16,76544 168,52444 397,27487 534,42041 62,76883 [8] 200,08041 56,51281

R specific loops

Ask for maximal value of elevation by county

Warning in FUN(X[[i]], ...): no non-missing arguments to max; returning -
Inf

[[1]] [1] 773,7128

[[2]] [1] 356,1018

[[3]] [1] 144,0656

[[4]] [1] 387,1294

Example: static map - locally stored data

First we need to prepare data

Inspect the data

CRS arguments: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

FIPS ISO2 ISO3 UN NAME 0 AC AG ATG 28 Antigua and Barbuda 1 AG DZ DZA 12 Algeria 2 AJ AZ AZE 31 Azerbaijan

Base graphic system, static map

Base graphic system, improvement of static map

Base graphics

Zoom, additional improvements

Additional improvements

Base graphic, different adjustments

[1] “Polygon” attr(,“package”) [1] “sp”

Formal class ‘Polygon’ [package “sp”] with 5 slots ..@ labpt : num [1:2] 17 44,5 ..@ area : num 25,2 ..@ hole : logi TRUE ..@ ringDir: int -1 ..@ coords : num [1:5, 1:2] 13,9 20,2 20,2 13,9 13,9 42,5 42,5 46,5 46,5 42,5

Different map, base graphics

  • basline data the same

Improved map

Leaflet example, selected polygons/counties

  • Oroslavje and Lobor

Example: data from MS Access

  • Not spatial DB ..

galeb

Create ‘spacetime’ class object - package ‘spacetime’

plotKML version 0.5-8 (2017-05-12)
URL: http://plotkml.r-forge.r-project.org/
KML file opened for writing...
Writing to KML...
Closing  galeb_52_3d.st_600.kml

Example: data from PostGres/PostGIS

Examples of analysis, points neighbours (sp)

Nearest neighbourgs, k=2; library(spdep)

Visualisation of neighbours

Analysis of point pattern

[1] “SpatialPoints” attr(,“package”) [1] “sp”

Planar point pattern: 300 points window: rectangle = [13,556343, 19,210248] x [42,53421, 46,50464] units

real-valued pixel image 128 x 128 pixel array (ny, nx) enclosing rectangle: [13,55634, 19,21025] x [42,53421, 46,50464] units dimensions of each pixel: 0,0442 x 0,031019 units Image is defined on the full rectangular grid Frame area = 22,4484489889368 square units Pixel values range = [13,36395, 13,36395] integral = 300 mean = 13,36395

Visualisation of result, export is necessary

Library geosphere

When data are not in Cartesian system

Visualising results

Call external software for some algorithm

  • Example calling SAGA GIS geoprocessor
    • Possible to have several versions of SAGA, specify in the call
  • QGIS
  • ESRI …

End, part 2!